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Geometric frustration leads to complex phases of matter with exotic properties. Antiferromag- 
nets on triangular lattices and square ice are two simple models of geometrical frustration. We map 
their highly degenerated ground-state phase spaces as discrete networks such that network anal- 
ysis tools can be introduced to phase-space studies. The resulting phase spaces establish a novel 
class of complex networks with Gaussian spectral densities. Although phase-space networks are 
heterogeneously connected, the systems are still ergodic except under periodic boundary conditions. 
We elucidate the boundary effects by mapping the two models as stacks of cubes and spheres in 
higher dimensions. Sphere stacking in various containers, i.e. square ice under various boundary 
conditions, reveals challenging combinatorial questions. This network approach can be generalized 
to phase spaces of some other complex systems. 



1. INTRODUCTION 

One challenge to understanding disordered solids is the complex geometry of their phase spaces, including the 
relative positions and interconnections between the different metastable states. Phase spaces are usually too large 
and complicated to be directly studied. For example, an A^-particle system typically has a vast abstract 67V-dimension 
phase space {3N for position, 3iV for velocity). Here, we propose that some simple models of disordered solids, such 
as geometrical frustrated spin models, provide an ideal platform for phase-space studies. Their phase spaces can be 
mapped as nontrivial complex networks, so that the recently developed large tool box of network analysis [Ij [21 [3] can 
be used to understand phase spaces. On the other hand, these phase spaces provide a new class of complex networks 
with novel topologies. 

When a system has competing interactions, there is no way to simultaneously satisfy all interactions, a situation 
known as frustration. Frustration widely exists in systems ranging from neural networks to disordered solids. Frus- 
tration can also arise in an ordered lattice solely from geometric incompatibility [A - For example, consider the three 
antiferromagnetic Ising spins on the triangle shown in Fig. [l]\. Once two of them are antiparallel to satisfy their 
antiferromagnetic interaction, there is no way that the third one can be antiparallel to both of the other two spins. 
Frustration leads to highly degenerated ground states and, subsequently, to complex materials with peculiar dynamics 
such as water ice [S] , spin ice [5] , frustrated magnets [5] , artificial frustrated systems [7] and soft frustrated materials 
i- 

In geometrical frustrated systems, spins on lattices have discrete degrees of freedom, such that their phase spaces 
are discrete and can be viewed as networks. A node in the network corresponds to a state of the system. Two 
nodes are connected by an edge (i.e. a link) if the system can directly evolve from one state to the other without 
passing through intermediate states. Edges are undirected because dynamic processes at the microscopic level are 
time reversible. The challenge is how to construct and analyze such large phase-space networks. For example, how 
do we identify whether or not two nodes are connected? 



2. ANTIFERROMAGNETS ON TRIANGULAR LATTICES. 

The first model we consider is antiferromagnetic Ising spins on a two-dimensional (2D) triangular lattice [9^. For 
a large system with periodic boundary conditions, it has ~ g0.3237V3pi„ degenerated ground states where Ngpin is the 
number of spins [9 . For example, configuration 3A in Fig. [ijC is one ground state in the hexagonal area. We refer to 
pairs of neighbouring spins in opposite states as satisfied bonds, i.e., they satisfy the antiferromagnetic interaction. 
Since one triangle has at most two satisfied bonds (see Fig. [iJV), the ground state should have 1/3 of its bonds 
frustrated and 2/3 of its bonds satisfied [9j. If we plot only satisfied bonds, a ground state can be mapped to a 
random lozenge tiling [10], see configuration 3A in Fig. [ip. A lozenge is a rhombus with 60° angles. By colouring 
lozenges with different orientations with different grey scales, the tiling can be viewed as a stack of 3D cubes, or as a 
simple cubic crystal surface projected in the [1,1,1] direction [10 , see Fig.[lfD. 

The ground state has a local zero-energy mode, as shown in Fig. [l|3: the central particle can flip without changing 
the energy since it has 3 up and 3 down neighbours. The system can evolve via a sequence of such single spin flips, 
even at zero temperature. We call such a local zero-energy mode the basic flip. Any configuration change can be 
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FIG. 1: A phase-space network of cube stacks. (A): Three antiferromagnetic spins on a triangle cannot simuhaneously 
satisfy aU their interactions. (B); The central spin has three up and three down neighbours, so that it can flip freely without 
energy change. Satisfied bonds can be viewed as cubes. The +/- free spin flip corresponds to adding/removing a cube. (C): The 
2x2x2 cube stacks are stable against gravity along the [1,1,1] direction. Stack configurations have one-to-one correspondence 
to Ising ground states under 'hexagon boundary condition', e.g., see configuration 3A. In the right 3A configuration, the black 
lines are satisfied bonds forming rhombuses and the blue lines are frustrated bonds. In total, there are 20 legal stacks, i.e., 20 
nodes in the phase-space network. The network is bipartite, i.e., consisting of alternating red (even number of cubes) and black 
(odd number of cubes) states. 



viewed as a sequence of such basic flips. Recently, we directly observed such flips in a colloidal monolayer [H]. In 
the language of cubes, a basic flip is equivalent to adding or removing a cube, see Fig. [Tj3. By continuing to add or 
remove one cube from the stack surface, we can access all possible stack conflgurations in the large box. Thus, the 
ground-state phase space is connected by this 'hexagonal boundary condition'. The corresponding cube stacking in a 
large box is equivalent to the boxed plane partition problem in combinatorics [11]. The total number of ways to stack 
unit cubes in an box is given by the MacMahon formula: [T^] 

TT i + j + k-l _ H^{L)H{3L) 

l<i,j,k<L ■' ^ ' 



27^ ^ 
16 



when L — > oo, (1) 



where the hyperfactorial function H{L) = HLo The first several N„{L = 2,3,4,5, •••) are 20, 980, 232848, 
267227532,- • • (see the number sequence A008793 in ref When L = 2, all 20 ground-state configurations in 

Fig. [Tp have the same minimum possible energy, i.e., 12 frustrated bonds in 12 rhombuses. 



3. NETWORK PROPERTIES. 

The 20-node phase-space network shown in Fig.[lp can be constructed based on the following two facts: (1) Spins 
have discrete degrees of freedom, such that the phase space is a discrete network; (2) Any configuration change can 
be decomposed to a sequence of basic flips. Consequently, we can define an edge between two nodes if the two states 
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FIG. 2: Connectivity distributions. Histograms of the connectivities of ground-state piiase-space networks for (A) anti- 
ferromagnets on triangular lattices, and (B, C, D) square ices under different boundary conditions. (A): L = 4 (circles) and 
1/ = 3 (squares) cube stacks. (B): Spheres stacks in L = 6 (circles) and L = 5 (squares) tetrahedra. (C): Spheres stacks in 
L = 4 (circles) and L — 3 (squares) octahedra. (D): Sphere stacks in L = 6 (circles), L = 5 (squares) and L — 4 (diamonds) 
containers shown in Fig. Insets: semi-log plots. The curves in the main plots and insets show the best Gaussian fits. 



differ by only one basic flip (i.e., one cube), such that the system can directly change from one node to the other 
without passing through intermediate nodes. 

Numerically, we can handle networks only up to L = 4 stacks with Nn = 232848 nodes; nevertheless many general 
properties have emerged from such small systems. Figure [2]A shows the connectivity (i.e. degree) distribution 
of cube-stack networks. The connectivity, ki, is the number of edges incident with the node i. The connectivities 
of various frustrated systems appear to have Gaussian distributions (see Fig. |2|. This behavior is similar to that of 
small- world networks [31 [T3] and Poisson random networks ^ i3j and different from that of scale-free networks (31 [TS] • 

Other network properties, such as the diameter and the cluster coefficient, can be readily derived from the cube 
stack picture. The shortest path length between two nodes is simply the number of different sites among all the 
sites. The largest distance, i.e. the diameter of the network, is between the 'vacant' and the 'full' states. Here, 
we define the vacant state as no cube (i.e., vacant sites) and the full state as no vacant site (i.e., cubes). 
The networks have small- world properties |31 |T3] in the sense that the diameter, L^, is almost logarithmically small 
compared with the network size, ~ g^^pin ^ ^ r^^ie network is bipartite (see the black and red circles in Fig. ^p) 
because a cube stack comes back to its initial configuration only by adding and removing the same number of cubes, 
i.e., an even number of basic flips. Consequently, the cluster coefficient [3 , which characterizes the density of triangles 
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in the network, is 0. 

Spectral analysis provides global measures of network properties. For an iV„-node network, the connectivity (or 
adjacent) matrix A is an Nn x Nn matrix with Aij = 1 if nodes i and j are connected, and zero otherwise. Since edges 
in phase-space networks are undirected, A is symmetric and all its eigenvalues, A^, are real. The spectral density of 
the network is the probability distribution of these Nn eigenvalues: p(A) — -j^ J2f=i ^i^ ~ ^i)- P{^)'^ 9th moment, 

Mq, is directly related to the network's topological feature. Dq — NnMq — the number of paths (or 

loops) that return back to the original node after q steps [5]. In a bipartite network, all closed paths have even steps 
so that all odd moments are zero. Consequently, the spectral density is symmetric and centered at zero. The ith 
node with ki neighbours has ki ways to return back after two steps; hence, the variance cr^ = M2 = ^i/^n = k, 
where k = 2Nedge/Nn is the mean connectivity. We rescale the spectral densities by fc^/^ to the unit variance (see Fig 
[3]). The rescaled spectral densities of different frustration models collapse onto the same Gaussian distribution. By 
counting Dq, we show that spectral densities are Gaussian at the infinite-sized limit (see Section I of Supplementary 
Information (SI)). This distinguishes phase spaces from other complex networks. For example, the spectral density 
of a random network is the semicircle in Fig. [3] The spectral densities have triangular distributions for scale-free 
networks and irregular distributions for small- world, modular hierarchical and many real-world networks 1171 . 

Spectral analysis can also detect the network's community (or modular) structures |18j if there are any. The 
algorithm in ref. |18| identifies some relatively highly connected subnetworks (i.e., communities). However, we still 
observe a number of edges between subnetworks such that the whole phase space has to be considered as fully ergodic. 
Our simulation shows that the system can easily travel through the whole phase-space network via basic flips and will 
not be trapped in a local community for a long time. 

4. POISSON PROCESSES AND EQUAL PROBABILITY IN PHASE SPACES. 

The fundamental assumption of statistical mechanics is that the dynamic trajectory of a system wanders through 
all its phase spaces and spends the same amount of time in each equally sized region of the phase space. However, 
this 'equal a priori probability postulate' (essentially the same as the 'ergordic hypothesis' [19]) is not necessarily 
true, as Einstein noted [H]. How the system moves from one configuration to the next depends on the details of the 
molecules' interactions (e.g. nearest-neighbor antiferromagnetic interactions here); these microscopic dynamics may 
make some configurations more likely than others. 

Network analysis provides an opportunity to study ergodicity. Unlike billiards with deterministic trajectory, we 
assume the spin flipping is due to the random thermal motion and not depends on history. Thus the dynamical 
evolution of the system can be viewed as a random walk on its phase-space network. It still interesting to ask whether 
this random walk can uniformly visit each node given the complex topology of the network. In another word, whether 
the system can visit each possible microstate configuration under the complex constraint of local nearest-neighbor 
interactions. 

Random walks on a network are rather chaotic, and nodes with higher connectivities will be visited more frequently. 
Thanks to the theorem in ref. |20j . the mean visiting frequency for node i is ki/N^dgej which only depends on local 
connectivity, ki, and does not depend on the global structure of the network. Here, A^edge is the total number of edges. 
This theorem is a direct consequence of the undirectedness of edges. Although highly connected nodes are visited 
more frequently (~ ki), interestingly, the equal-probability postulate does not break down because the average time 
stayed at node i is ~ l/ki. Basic fiips are random and independent of history, meaning that it is a Poisson process. We 
define the fiipping probability of a basic fiip within a unit of time as which is the intensity of the Poisson process. 
In Poisson processes, the time interval between fiips (i.e., the staying time) has an exponential distribution, e~'^*, and 
the mean staying time is Note that multiple fiips will not fiip exactly simultaneously because time is continuous. 
Therefore we do not need to worry about possible illegal configurations if neighbor free spins fiip simultaneously. For 
a node with connectivity k, the superposition of k Poisson processes is still a Poisson process with intensity kiy and, 
consequently, the mean staying time is l/(kiy). A random walker has higher frequency (~ k) to visit a high-fc node, but 
will stay there for a shorter time (~ 1/^), so that the equal-probability postulate is recovered. Boltzmann assumed 
that molecules shift from one microscopic configuration to the next in such a way that every possible arrangement 
is equally hkely, i.e., all edges have the same weight. We find that the equal-probability postulate still holds if edges 
have different weights (see Appendix B), which, for example, can represent different potential barriers in complex 
energy landscapes in phase spaces. 
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FIG. 3: Spectral densities of phase-space networlcs. Variances are rescaled to 1 by A' = A/fc 2 . Black curve: Gaussian 
distribution ''^/%/27r. Dasiied curve: Wigner's semicircle law for random networks. p(A) = ^/Aa'^ — A^/(27r(T^) if |A| < 2a 
and zero otherwise. The variance is also rescaled to 1. Red curve: the spectral density of the 980-node network of L = 3 cube 
stacks. Blue curve: the 7436-node network of L = 6 sphere stacks in a tetrahedron, i.e., 7x7 square ice under the domain wall 
boundary condition. Green curve: 7782-node network of 2 x 5 square ice under the free boundary condition. Their Gaussian 
fits are indistinguishable from the black curve. 



5. SQUARE ICE. 

We further study another frustration model called square ice to identify the more general properties of phase spaces. 
Square ice is the two-dimensional version of water ice as shown in Fig. |9] It can be viewed as jigsaw tiling [22^ or spin 
ice [6j III |22l |23] (see Figs. |4]A_,D). Oxygen atoms are represented by vertices and the relative directions of hydrogen 
atoms are represented by arrows. The ground state of the system follows the ice rule, i.e., each vertex has two incoming 
and two outgoing arrows. It is also known as the six-vertex model since each vertex has six possible configurations 
(i.e., six types of jigsaw tiles). For a vertex associated with four ferromagnetic spins, frustration is inevitable (see the 
example in Fig. |4p). 

Flipping a closed loop of arrows from clockwise to counterclockwise (or vice versa) does not break the ice rule. The 
smallest four-spin loops in Figs.|4|D,G are labeled in red (clockwise) and yellow (counterclockwise). They are basic 
flips since any configuration change can be decomposed as a sequence of such flips |24]. Similar to cube stacking, 
all the legal configurations of square ice are connected via basic flips [21] . Consequently, the phase-space network of 
square ice can be constructed. 
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The square ices in Figs. |4j\,D have domain wall boundary conditions (DWB) as shown by the black arrows in 
Fig. |4|D. There is a one-to-one correspondance between jigsaw tiling with DWB and alternating sign matrices (ASM) 
[22] (see Fig.|4j\). ASM are square matrices with entries or ±1 such that each row and column has an alternating 
sequence of +1 and -1 (zeros excluded) starting and ending with +1. The number oi n x n ASM is [22] 

27 \ ^ 

— I when n oo, (2) 

ley 

i.e., the number of nodes of the phase-space network of an n x n square ice with DWB. 



6. MAPPING SQUARE ICES TO SPHERE STACKS. 



Mapping 2D triangular antiferromagnets to 3D cube stacks greatly simplifies the picture of the phase space and 
allows combinatorial analysis to generate quantitative results such as Eq. [l]and Gaussian spectral densities. Here, we 
show that square ices can be mapped to 3D close-packed spheres in face-centered cubic (FCC) lattices. Each square 
plaquette in Figs. |4]3,G is assigned a height [531 US] based on the rule shown in Fig. [4j3: When walking from the 
plaquette with height h to its neighbor, the height increases by 1 if it crosses a left arrow and decreases by 1 if it 
crosses a right arrow. The ice rule guarantees that the height change around a vertex is zero and h is independent 
of the path along which it was computed. From the minimum and maximum possible heights, we found that DWB 
yields a stack of building blocks in a tetrahedron (see Fig. S5 of SI). A plaquette can be flipped only when its four 
neighbor plaquettes have the same height. Since each building block is 'supported' by four underneath blocks in 
an effective 'gravity field', the stack can be viewed as an FCC lattice along the [100] direction, see section III of 
SI. Thus, the stacking blocks should be rhombic dodecahedra, which are primitive unit cells of FCC lattices. FCC 
lattices can be conveniently represented by close-packing of spheres. Sphere stacks in side length L tetrahedra have 
one-to-one correspondence to (L -I- 1) x {L + 1) square ices with DWB (see Fig. S5 of SI). The stack of red spheres 
in Fig. corresponds to the configurations in Figs. |4j\,D. The physical heights of the red spheres on the top surface 
are the heights of the corresponding plaquettes in the square ice. At the interface between the red spheres and the 
yellow vacant sites, the four removable red spheres on the top surface correspond to the red plaquettes and the four 
addable yellow sites correspond to the yellow plaquettes in Fig. [4p. Similar to the cube stacking, here, a basic fiip 
from counterclockwise to clockwise (or vice versa) corresponds to adding (or removing) a sphere. By adding spheres 
from the vacant state shown in Figs. [TO|\,G, we can generate all possible stack configurations, i.e., all the nodes of 
the phase-space network. Similar to the cube stack case, we can construct the phase-space network of sphere stacks 
by adding an edge between two nodes if the two stacks are different by one sphere. 



7. NETWORK PROPERTIES OF SPHERE STACKS. 



We numerically studied phase-space networks of small square ices under various boundary conditions. Our largest 
network contains 2068146 nodes and 13640060 edges (4x5 ice under free boundary conditions). All networks have 
the small-world property. Their connectivity distributions in Figs. [2}3,C,D and the spectral densities in Fig. [3] are 
similar to those of cube stacks. Apparently, both cube-stack and sphere-stack phase spaces have Poisson processes 
with equal probability and Gaussian spectral densities. 



8. BOUNDARY EFFECTS. 



Stacks in higher dimensions provide a vivid means for qualitative visualization of the boundary effect, which has not 
been well understood in geometrical frustration [26J . One peculiar property of geometrical frustration is that boundary 
effects often percolate through the entire system even in the infinite-sized limit [ 251 HZ] . This can be visualized from a 
typical sphere stack in the i = 100 tetrahedron shown in Fig.j^Ji^, which has a central disordered region and four frozen 
(ordered) corners known as the arctic circle phenomenon 28J. The disordered region is not uniformly random since 
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FIG. 4: Square ice. (A): A 5 x 5 square ice under the domain wall boundary condition. Each jigsaw tile can be viewed as a 
water molecule with one oxygen atom in the center and two hydrogen atoms at the two bulges (see Fig. S4 of SI). By assigning 
vertical tiles to be 1, horizontal tiles to be -1 and the other four types to be 0, a 5 x 5 alternating sign matrix |22] is obtained. 
(B): The height rule used in (D) and (G). (C): Four magnets placed at a cross inevitably have frustrations. (D): The spin ice 
mapped from (A). The arrows represent bulge directions in (A). The blue arrows may flip under the ice rule. Each plaquette 
is assigned a height based on the rule in (B). The upper left corner is defined as height zero. Basic flips (i.e., four-arrow loops) 
are labeled in red (clockwise) and yellow (counterclockwise). (E): The corresponding sphere stack of (D). Yellow spheres are 
vacant sites. (F): A typical sphere stack in an L = 100 tetrahedron. The sphere centers are connected so that it appears to be 
a stack of polyhedra. (G): A spin ice configuration in an Aztec diamond area under the constant- height boundary condition. 
(H): The corresponding sphere stack of (G) in an octahedron. 



different positions have diflFerent mean surface curvatures and entropy densities [57] (see Appendix C). Consequently, 
the infinitely large limit under DWB cannot be called the thermodynamic limit due to the lack of homogeneity. 

Different boundary conditions in square ice correspond to different container shapes in sphere stacking. For example, 
the boundary condition shown in Fig. |4]GI corresponds to sphere stacks in an octahedron (see Fig. |4j3) because the 
lowest possible heights form an inverted pyramid (i.e. the container) and the highest heights form an upright pyramid 
(i.e. the lid). Appendix D shows another boundary condition whose container and lid have different shapes. At 
a given boundary condition, the lid and the container form an interesting pair of dual surfaces. Some boundary 
conditions do not have the arctic circle phenomenon, as illustrated by the sphere stacking in Appendix C. 

We found that phase spaces are ergodic under free or fixed boundary conditions, but nonergrodic under periodic 
boundary conditions whose networks consist of disconnected subnetworks. As an example. Fig. [5] is the phase-space 
network of the 2x3 square ice wrapped on toroid, i.e., under periodic boundary conditions. It contains two nontrivial 



(12-node) subnetworks and 20 trivial isolated nodes. The corresponding 44 configurations are shown in Fig. 14 For 
TO X n periodic square ice wrapped on a toroid, we show that its phase space contains 2"+^ + 2"'+^ — 4 trivial isolated 
nodes, (m — 1) x (n — 1) nontrivial subnetworks and the smallest non-trivial subnetwork has nodes (see 

Appendix E). These results are confirmed numerically. 
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FIG. 5: The phase-space network of 2 x 3 square ice under periodic boundary conditions. There are 44 possible 
states (see their detailed configurations in Figs. S9A,B of SI). States 1 to 12 are connected via basic flips; states 13 to 24 are 
connected; and states 25 to 44 do not contain any basic flips, thus they are isolated nodes. 

9. DISCUSSION AND OUTLOOK. 

We build novel connections between geometrical frustration, combinatorics (e.g., plane partition and sphere stack- 
ing) and complex networks to exploit open questions and analysis tools from these fields. Other frustration models, 
such as triangular and kagome ices, antiferromagnets in 2D kagome and 3D pyrochlore lattices [5111 131j . have height 
functions and basic flips so that their phase-space networks can be similarly constructed. In principle, these models 
can be mapped to polyhedra stacking in higher dimensions, so that their rich symmetries and boundary effects become 
more transparent. 

Quasicrystals can also be mapped to higher-dimensional lattices. Projecting the high-dimensional lattices to lower 
dimensions could result in periodic lattices (i.e., geometrical frustration) at certain projection angles, or aperiodic 
structures (i.e., quasicrystals) at other angles. Phasons in quasicrystals correspond to basic flips in geometrical 
frustration [27], thus similar phase-space analysis may be applied to quasicrystals. In fact, the infinitely degenerated 
where N —> oo) ground states in both geometrical frustration and quasicrystals are essentially metastable 
states since the third law of thermodynamics dictates that the true ground state of real materials must have finite 
degeneracy. Network analysis may provide a possible approach to understanding the observed glassy dynamics in 
frustrated systems [S]. At finite temperatures, phase-space networks can be similarly constructed. The nodes are all 
configurations on the hypersurface in the phase space determined by the conservation laws. Configurations change 
via basic flips and diffusion of thermal excitations [8l[32]. These motions are represented by edges. The weight of 
each edge can be assigned by a Boltzmann factor or defined by the physical details of the real system [35] . Height 
representation can be recovered by assigning vector heights [33], so that systems at finite temperatures might be 
mapped to stacks in even higher dimensions. 

Compared with intensively studied social networks, information networks, biological networks and technological 
networks [2], phase-space networks belong to a new class with unique Gaussian spectral densities. A large tool box ^ 
has been developed in the recent decade to study network dynamics, correlations, centrality, community structures, 
fractal properties [Mj, coarse graining [3S], etc. These tools can be readily applied to phase-space studies. In 
particular, phase spaces might have fractal structures because stacks of cubes or spheres have self-repeating patterns 
on various length scales. This may cast new light on the highly controversial Tsallis's nonextensive entropy [36l I37j . 
which is based on the assumption that nonequilibrium systems have fractal phase spaces. To date, a real example to 
support this assumption has not been available. Indeed, geometrical frustrated ground states share the same features 
as the long-range interacting systems typically discussed in the context of Tsallis entropy. One example is boundary 
effects percolating through the entire system so that the system is not uniform at the infinite-sized limit and cannot 
be viewed as a simple sum of its subsystems (i.e., non-extensive). 

In statistical physics, the two models we studied here are considered as exactly solvable [38] under periodic boundary 
conditions at the infinite-sized limit. Combinatoric analysis, although challenging, provides an alternative approach to 
yield exact results about finite systems and at various boundary conditions. Cube stacking (i.e., rhombus tiling or plane 
partition [22]), naturally appears in many chemical and physical problems, such as counting benzenoid hydrocarbons, 
percolation, crystal melting and string theory |45j . In contrast to the intensively studied cube stacking, sphere stacking 
has not been explored. Only some combinatoric properties of sphere stacking in tetrahedra are available since we 
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can map them to ASM. Our numerical calculations show that there are 2, 7, 42, 429, 7436, 218348 ■• • ways to pack 
spheres in L = 1, 2, 3, 4, 5, 6, • • • tetrahedra; and 2, 18, 868, 230274, • • • ways in L = 1, 2, 3, 4, • • • octahedra. The former 
number sequence (i.e., sequence A005130 in ref.fl^) is given by Eq. [2] while the formula for the latter is not available. 
Moreover, many questions studied in cube stacking can be asked about FCC sphere stacking. For example, how many 
ways are there to pack N spheres into a tetrahedron? Is there a similar generating function as cube stacking for 
sphere stacking in a tetrahedron [T5]? What is the ensemble-averaged surface in Fig. |6jV, i.e., what are the entropy 
density distributions at the infinite-sized limit [271 HD]? These questions can also be asked about other container 
shapes. Furthermore, square ice has one-to-one mappings to other 2D models, such as three-color graphs, dimers, 
fully packed loops, etc. [23]. It also has one-to- multiple mapping to the domino tiUng [IT]. Sphere stacking provides 
a simple 3D picture and casts new light on these 2D models. 
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APPENDIX A: PROOF OF THE GAUSSIAN SPECTRAL DENSITY 



The characteristic function, i.e., the Fourier transform of the probability function, uniquely describes a statistical 
distribution. It can be written as a series of moments of the distribution. Hence, to prove that the spectral density 
is Gaussian, we only need to show that all orders of the moments are the same as those of a Gaussian distribution. 
For a Gaussian distribution centered at 0, its odd moments are zero and its even moments (of order q) are M{q) — 
2i/^q/2)\ ^'^ — il ~ l)-'"''; whcrc (T^ is the variance. For an A'n-node undirected network, D{q) = NnM{q) is the 
number of directed paths that return to their starting node after q steps [3j. We count Dq by stacking cubes/spheres 
and show that M{q) = D{q)/Nn follows the Gaussian M{q). 

The phase-space networks are bipartite since walking an odd number of steps (i.e., adding/removing cubes/spheres 
an odd number of times) cannot return back to the original state. Consequently, all odd moments are zero, i.e., the 
distribution is symmetric and centered at 0. D{2) is the number of ways to have one basic flip, /i, and its reverse flip, 
/i. Given a stack configuration, i, with ki available basic fiips, i.e., node i with connectivity ki in the phase-space 
network, its D{2)i = ki. Thus, the total D{2) = J^i^^i = ^nk where k = 2Nedge/Nn is the mean connectivity. 
Compared with the second moment, M(2) = cr^, we have k — a^. D{A) is the number of ways to have two basic fiips, 
/i, /2, and reverse fiips, /i, f2- Subscripts 1 and 2 denote the time order. Given /i and /2, typically there are three 
ways to arrange them in legal order: /1-/1-/2-/2, /1-/2-/1-/2 and fi-h-h-fi- Note that the reverse fiip, Jj, must 
be later than fj. fj represents either adding or removing a cube/sphere. If /i and /2 are fiips of the same spin or 
neighbor spins, they are not independent so that only /1-/1-/2-/2 is legal. However, the probability of such a case 
approaches in infinitely large systems such that we can neglect the 'interference' between basic fiips in large systems 
and assume that all fiips are independent. Next, we consider how many choices of /i and /2 we have. Given the initial 
state, i, /i has ki choices and /2 has kn choices. Here, kn is the connectivity of a node after walking one step away 
from state i. In large systems, a dominant number of states have large k diverging with the rough surface area, ~ L^. 
Hence, kn ~ ki. Moreover, the dominant number of states is close to the mean surface, such as Fig. [6]A_ under the 
domain wall boundary condition. The surface shape distribution peaks around this maximum possible surface and 
becomes like a Dirac delta distribution when approaching the infinite-sized limit |27j . The probability distribution of 
normalized connectivity approaches a Dirac delta distribution as well (see Fig. |7 and its caption). Thus the leading 
term in kiku is k^. At the infinite-sized limit, there are kiku ~ fc^ choices of /i and /2, and three ways to fiip them 
in different time orders; therefore, -D(4) = 3fc^. 

Similarly, we count D{2n) by considering 2rt fiips, /i, /i, /2, /2, • ' ' i In, fn- They are placed in a 2?T,-long sequence 
in time order. First, /i must be placed at step 1. Then, there are 2n — 1 choices for placing /i. Then, /2 must be 
placed at the earliest available step (i.e., step 2 if /i is not occupying that step). Then, /2 has 2n — 3 choices. Thus, in 
total, there are (2n— 1)!! legal sequences. Note that this is accurate because a finite number of /j's are diluted enough 
to be considered as independent in an infinitely large system. Next, we consider how many choices of fj's there are. 
Given the initial state, i, /i has ki choices, /2 has kn choices, ... /„ has choices. Here, kij is the connectivity 

of a node after walking j steps away from the initial node, i. kij depends on the pathway of the j steps and is not a 
constant. In total, there are 0^=0 % choices for {/i, /2, • • • , /«} if node i is chosen as the starting point. When the 
system size L 00, there will be Nn{l — S) states {S — > 0) with large connectivity, ki ^ L^. In cube stacks, adding one 
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FIG. 6: (A): The mean sphere stack surface in an L = 100 tetrahedron. The surface is obtained by averaging over 10^ stacks at 
equilibrium. A typical surface at equilibrium is shown in Fig. 4F of the main text. (B): The corresponding probability of basic 
flips measured from 10^ step simulation. The probability distribution appears to be a hemisphere. (C): The flipping probability 
in (B) represented by brightness. The bright non-frozen region is circular, which agrees with the arctic circle theorem [2S] . 




k' 



FIG. 7: The probability distribution of normalized connectivity, k' = k/kmax, where the maximum connectivity is kmax oc 
Nspin OC L^. This figure is normalized from Fig. 2 in the main text. Black curves: cube stacks in L = 4 (thick curve) and 
1/ = 3 (thin) boxes. Red curves: sphere stacks in L — 7 (thick) and L — 6 (thin) tetrahedra. Blue curves: 2D square stacks 
(i.e., 2D sphere stack as shown in Fig. 15 \) with L = 11, 10, 9, 8 (thicker curves for larger L). The peaks are in the middle bin 
at 0.5, i.e., k = kmax/2 has the highest number of stack configurations. In the infinite-sized limit, the normalized distributions 
will be asymptotic to a specific functional form. For a histogram of Fig. 2 in the main text before the normalization, the peak 
height, H, increases exponentially with the system size, while the heights at kmax are always 1 or 2 (for example, see Fig. [sf. 
Thus, their ratio 2/H — > in infinitely large systems. This indicates that the asymptotic distribution is indeed a Dirac delta 
function as described in ref . [37] . 



cube can change A: by 3 at most because one cube is supported by three underlying cubes. If the shortest path between 
two nodes has j steps, their connectivity difference is \Skij\ < 3j. In sphere stacks, one sphere is supported by four 
underlying spheres; thus, \5kij\ < 4j. Therefore, when L oo, j ^ L and 01=0 ~ 

by dropping the high-order terms. The last step uses the fact that the surface shape distribution becomes a delta 
distribution when approaching the infinite-sized limit. 

Combining the above results, D{2n) ~ X^ilTiC^'^ — 1) •■11^=0 — ~ ^)-^nk", which becomes exact at the 
infinite-sized limit. Since D{2) — Nnk — NnCr"^ , the 2nth moment of the eigenvalue distribution M{2n) — D{2n)/Nn — 
{2n — 1)!!(T^" is identical to the 2nth moment of a Gaussian distribution. Therefore, spectral densities of phase-space 
networks are Gaussian at the infinite-sized limit. In fact. Fig. 3 in the main text shows that spectral densities are 
already very close to the Gaussian distribution when systems are small (~ 10'^ nodes). 
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FIG. 8: The two highest connectivity (k = 13) states in L 
in L = 4 cube stacks. 



3 (odd) cube stacks and the highest connectivity (fc — 24) state 



APPENDIX B: DYNAMICS IN PHASE-SPACE NETWORKS WITH WEIGHTED EDGES 

In the main text, we show that trajectories spend equal amounts of time at every node on average. This can be 
easily generalized to networks with weighted edges, which, for example, can represent different potential barriers at 
finite temperatures. We replace an edge, i, with weight Wi with Wi equivalent edges. By applying such replacement to 
all edges, we get a new network with all equally weighted edges. The new connectivity for node i is = Wij. 
For the same reason as shown in the main text, the mean staying time at node i is oc 1/fc™. The visiting probability 
is oc k'!f by generalizing the theorem in rcf.[20. to weighted networks 42]. In fact, the proof in ref. ^20^ can be directly 
applied to weighted networks. For a weighted network, the connectivity matrix, Aij = Wij, where Wij is the weight of 
the edge between nodes i and j. The weighted connectivity is k^" — Aij. The transition probability from node i 
to node j is Aij/kf . Suppose a walker starts at node i at time t — Q. Then, the master equation for the probability, 
Pij, to find the walker at node j at time t is given by [20] : 



3t-l 



Jt-1 



The transition probability, Pij (t) , from node i to node j in t steps can be explicitly expressed by iterating Eq. [s] 

p -m- V An jVz2 ...^V^ 

ji,---,jt-i * -Ji •''-i 
By comparing the expressions of Pi^j{t) and Pj^i{t), we get 

krp,^,{t) = kjp.^s)- (5) 

We define the stationary probability, as t ^ oo. Eq. |5] implies that kiP°° — kjP°° and, consequently, we obtain 

uw \^ A ■ . 

poo _ l^i __ l^i ^'ij ff.\ 

A random walker visits node i at frequency P°° ^ kf and stays at node i for the time period ~ 1/^i" on average. 
Thus the walker spends the same amount of time at each node. 



APPENDIX C: SQUARE ICE AND SPHERE STACKS 

Figure [9] shows a 5 x 5 square ice with the domain wall boundary condition (DWB) corresponding to Figs. 4A, 
D, E. in the main text. Figure 10 shows 6x6 spin ices with DWB and their corresponding FCC sphere stacks in 
an L = 5 tetrahedron. Flipping all counterclockwise four-spin loops to clockwise corresponds to adding one layer of 
spheres. The tetrahedron emerges from the maximum packing, i.e., the full state. 

A square plaquette can flip only when its four neighbors have the same height, which indicates that one building 
block in the 3D stack should be supported by four building blocks underneath. Thus, the stack can also be viewed as 
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FIG. 9: (A) Square ice with a domain wall boundary condition. Water molecules are frozen into a square lattice. This 
configuration corresponds to Figs. 4A, D, E in the main text. 

a body-centered cubic (BCC) lattice ^5] . Both FCC and BCC stacking have the same combinatoric properties since 
they are only different by a stretch (see a 2D analogy in Fig. [TsjA). FCC is better than BCC because (1) FCC can 
be viewed as close-packed spheres in simple container shapes; (2) FCC has simple polyhedra stacking under gravity. 
The Wigner-Seitz cell of an FCC lattice is a rhombic dodecahedron, which is supported by four rhombic dodecahedra 
underneath, while the Wigner-Seitz cell of a BCC lattice can be supported by one block underneath since it has a flat 
square on the top. 

APPENDIX D: BOUNDARY EFFECTS 

We use the 'alternating boundary condition' shown in[TT|A, B to elucidate the lid-container duality. Given a bound- 
ary condition, the minimum (or maximum) possible heights can be directly written out, for example see Figs. |11[ \, 
B. These heights define the lid and the container. The lid and the container have different shapes. The container 
contains multiple height minima and the lid contains one height maximum (see the colored squares in Figs. [Tl] ^, B). 
They form an interesting pair of dual surfaces: using one as the container, the other will emerge as the surface of 
the highest 'sand pile' of small spheres. In fact, given a fixed boundary condition, the container and the lid are dual 
surfaces because, if we reverse the height rule in Fig. 4(B), the container and the lid switch roles. Packing spheres in 
the container is equivalent to packing buoyant spheres in the corresponding lid. The height difference between a lid 



and container is a pyramid; for example, see Figs. 11 and 13 



Some boundary conditions do not exhibit the arctic circle phenomenon shown in Fig. [6j Their disordered region 
may not have a circular shape or may not even have a frozen area [24j. The container shape provides an intuitive 
understanding about how boundaries affect the disordered region. In a tetrahedron, the largest horizontal cross 
section is a square in the middle height that is circumscribed by the disordered region. In an octahedron, the largest 
horizontal cross section is the total square ice area so that there is no frozen region under the boundary condition 
shown in Fig. 4G of the main text. This is confirmed by our simulation. Other boundary conditions may lead to 
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FIG. 10: The one-to-one mapping between square ice and sphere stacks. (A-F): 6x6 square ices. The domain waU boundary 
condition is labeled in grey. The height of each square plaquette is labeled with a number ranging from to 6. The upper 
left corner is defined as zero height. Other heights are generated by the height rule in Fig. 5B in the main text. (A): the 
vacant state with the lowest possible heights. Heights in (F) define a container consisting of two yellow triangles shown in (G). 
Counterclockwise four-spin loops are labeled in yellow. Flipping all of them in (A) leads to the configuration in (B). Clockwise 
four-spin loops are labeled in red. Flipping the yellow plaquettes results a series of states shown in (C, D, E, and F). (F) has 
the highest possible heights without yellow plaquettes. These heights define a lid, which is an upside-down container in (G). 
All legal square ice configurations can be generated by flipping yellow plaquettes from the vacant state or, reversely, by flipping 
red plaquettes from the full state. A plaquette can be flipped only when its four neighbour plaquettes have the same height. 
The height of each plaquette is the physical height of the corresponding spheres on the top surface of the stack. Each basic flip 
can be viewed as adding or removing a sphere. (H-L): the red sphere stacks corresponding to (B-F). The yellow spheres are 
vacant sites. 



the non-circular disordered region. For example, the flower shape observed in ref. |24j (see Fig. 12 3) is a direct 



consequence of the container shown in Fig. 12 \_ 



The ensemble average over random stacks results a mean surface shown in Fig. SIA. At the infinite-sized limit, 
dominate states are close to this mean surface [57], i.e., the typical surface in Fig. 4F approaches the mean surface 
in Fig. SIA when L — > oo. In the space of the stack surface, the distribution peaks around this maximum possible 
surface and becomes more and more like a Dirac delta distribution when approaching the infinite-sized limit [27 . The 
local gradient of the surface determines the density of the basic fiips, i.e., the density of configurational entropy so 
[27j . In Fig. SIA , the sq is zero in the frozen areas and continuously varies to reach its maximum value near the 
center of the square ice, with a non-zero gradient everywhere except near the center [27|. Consequently, the infinitely 
large limit of the DWB cannot be called the thermodynamic limit due to the lack of homogeneity. In contrast, 
the boundary condition in Fig. 4G has the thermodynamic limit since the limiting surface in the octahedron is flat 
everywhere. The flat surface has the maximum possible sq (for example, see Fig. S3) so that its spatial averaged. 
So, is as high as that of the free boundary condition. The boundary condition in Fig. 4G is a subset of the periodic 
boundary condition, so that the periodic boundary condition has the same sq ^ the free boundary condition at the 
infinite-sized limit. This explains why the sg calculated from the periodic boundary condition [43^ agrees so well with 
the experimental results on water ice obtained under the free boundary condition [33] . When the height difference of 
a fixed boundary is comparable to L, the limiting surface is not fiat and sq is smaller. For example, the zero-point 
entropy of L x L square ice in the infinite-sized limit under DWB is So = ksL'^ In Nn = fcBln(^27/16) [41 based 
on Eq. 2 in the main text, which is smaller than kg L'^ ln( -^64/ 27) under the periodic boundary condition [29 . For 
cube stacks, sq = ks ln(-\/27/16) under the hexagonal boundary condition based on Eq. 1 in the main text, which is 
smaller than 0.323fcB under the periodic boundary condition [9 . 

The typical stack configuration in a tetrahedron or octahedron is ^ 50% filled because the lid and the container have 
the same shape. When they have different shapes, a typical stack configuration may not be ~ 50% filled. Figure 12 \ 
is the averaged surface in an L = 100 container. The total volume has L{L^ — l)/6 spheres, and ^42% of the volume 
is filled with spheres on average. Note that the largest horizontal cross-section is at = a/2/3 with a corresponding 
filled fraction of 4/9. 
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FIG. 11: 6x6 square ice with the 'alternating boundary condition'. (A): The maximum possible heights, i.e., the lid, contains 
one clockwise basic flip labeled in red. (B): the minimum possible height, i.e., the container, contains multiple counterclockwise 
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5=. 20 40 60 SO 100 



FIG. 12: (A): The ensemble-averaged height surface for L = 100 square ice with alternating boundary condition shown in 
Fig. [11]^, B. The heights are rounded off to integers to show equal-height contours. The container shape is shown in Fig. |llp . 
(B): The flipping probability represented by the brightness. 



APPENDIX E: SQUARE ICE WITH PERIODIC BOUNDARY CONDITION 



Here, we use the 2x3 square ices to illustrate that periodic boundary conditions result in disconnected phase-space 



networks. Figure 14 shows the 44 configurations of the 2x3 square ice wrapped on a toroid. Note that this periodic 
boundary condition is for spins, not for heights. The upper left corner is defined as zero height. The 12 configurations 
in Fig. [l"4|V are connected by basic fiips and form a 12-node cluster as shown in Fig. [Tip. The other 12 configurations 
in Fig. [14) 3 form another 12-node cluster in Fig. [T4p. The height difference between the top corners and bottom 
corners is -1-1 in Fig. 14 \ and -1 in Fig. [14)3. Note that the four corners are essentially the same plaquette on the 
toroid, so they must be either all inside or all outside of a loop, such as the one shown in Fig. [Tip. Consequently, 
the height difference between the corners cannot be changed by flipping a closed spin loop. Therefore, the nodes in 
Fig.[T4?\_ and B form two disconnected clusters. The 20 isolated nodes in Fig. [Tip correspond to the 20 configurations 
in Fig. Tip, which contain no basic fiips. 

We generalize the above results to the m x n periodic lattice and prove that it has (m — 1) x (n — 1) non-trivial 
clusters and 2"+^ -I- 2™+^ — 4 isolated nodes. Unlike fixed boundary conditions, after walking along a closed loop 
in the x or y direction on a toroid and coming back to the original plaquette, the height may change. Such height 
differences, Ah^ and Ahy, uniquely characterize each disconnected subnetwork. Consider an arbitrary plaquette on 
a toroid. We unwrap the lattice onto a plane so that this plaquette is at the upper left corner with the height defined 
as 0, e.g., see Fig. [Tip. As shown in Fig. [Tip, any zero energy fiip of a spin loop cannot change the height differences 
between the four corners since they are essentially the same plaquette on the toroid. Consequently, configurations 
with different Ah^ or Ahy cannot be connected by basic flips. On the other hand, if configurations have the same 
Ahx and Ahy, they must be connected because they have essentially the same fixed boundary condition after being 
unwrapped onto a plane (see Fig. 14 Here, we use the fact that all legal configurations at a fixed boundary condition 
are connected by basic fiips [Hj. Since they are connected, we can choose the configuration whose bulk spins are 
along the boundary spins to represent each subnetwork (see examples in Figs. lAJ, F). Next, we consider the number 



of representative configurations, i.e., the number of subnetworks. If a representative configuration has no basic fiips, 
all of its horizontal spins or vertical spins have to be along the same direction as shown in Fig. [Tip. The 'energy 
barriers' between these states are large for large systems because m or n spins need to fiipped simultaneously in order 
to change from one state to another without breaking the ice rule. If all horizontal spins are leftwards (or rightwards), 
there are 2™ configurations for vertical spins (see Fig. lA'J). If all vertical spins are upwards (or downwards), there 
are 2" configurations for horizontal spins. In total, the four configurations are double counted so that there are 
2"+! _|_ 2™+i — 4 isolated nodes, i.e., configurations without basic fiips. Next, we consider nontrivial clusters with 
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FIG. 13: The lid, container and their height difference of a tetrahedron. 



multiple nodes. The corner height difference, lS.hx, has m — 1 possible values, and A/iy has n — 1 possible values 
(see Fig. 14?"), so that there are (m — 1) x (n — 1) nontrivial subnetworks in total. The representative configurations 



of the six subnetworks shown in Fig. are chosen to have the lowest possible heights, i.e., vacant containers for 
sphere stacking. Each configuration is characterized by one basic flip labeled as yellow squares, i.e., the lowest point 
of the vacant container. Apparently, there are (m — 1) x (n — 1) positions for a yellow square, i.e., (m — 1) x (n — 1) 
subnetworks. This result confirms that the zero-point entropy sq of the whole network is the same as that of the largest 
subnetwork under the constant-height boundary condition (see section IV) because (m— 1) x (n — 1) is logarithmically 
small compared with the total number of configuration ~ g^spi" ^ e™ ". 



Next, we show that the smallest nontrivial cluster has 



(m+K-l)! 



nodes. In Fig. 



14 ^, the two middle configurations 
subnetworks. When the yellow 



(n-l)!(m-l)! 

represent 132-node subnetworks and the other four configurations represent 60-no'3e' 
square is at the corner, the height function indicates that the container shape is a tilted 2D container rather than a 
3D container. Thus, the number of spheres packing in such a container is much smaller than that in 3D containers 
whose lowest point (the yellow plaquette) is not at the corner. To count the number of states in a tilted 2D container, 
we first consider the simple case in Fig. [TllA. C onfiguration 1 in Fig. 14\_ is the representative state with the lowest 
height of -3. Configurations 1 to 4 in Fig. |14|\ have the same boundary spins so that we can view them as the 2D 
sphere stacks in the same 1 x 3-sized 2D rectangle. Such a blue 1x3 container has three possible positions relative 
to the zero height plaquette (see configurations 1, 5 and 9 in Fig. \Ak). In total the subnetwork has 3 x 4 = 12 
nodes. It is easy to generalize this counting to m x n square ice on a toroid. There are m possible positions for 
the m X {n — l)-sized rectangle. With 2D sphere stacking in an a x 6 container, there are C^_|_^ — {a + b)l/{albl) 



configurations (see Fig. 
nontrivial subnetworks 



15 



and its caption). Consequently, there are mC^ 



-1)! 



nodes in the smallest 



m+n-l („_i)!(m-l) 

We confirm the above results numerically. Our numerical results also confirm the number 



sequence A054759 in ref.[T3] for the n x n square ice under periodic boundary conditions. 



[1] Albert, R. & Barabasi, A.-L. Statistical mechanics of complex networks. Rev. Mod. Phys. 74, 47-97 (2002). 
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FIG. 14: Configurations and tlie pliase-space network of 2 x 3 square ice wrapped on a toroid, i.e., witfi tlie periodic boundary 
condition. (A): Configurations 1 to 12. The upper left plaquette is defined as zero heiglit and other heights are derived from 
it by using the height rule in Fig. 4B in the main text. (B); Configurations 13 to 24. (C): 20 configurations that contain no 
basic flip. They are categorized into four types: all horizontal spins are (1) leftwards; (2) rightwards; all vertical spins are (3) 
upwards; (4) downwards. Four configurations are double counted, so the total number of configurations is 2 x 2^ -I- 2 x 2'^ — 4. 

(D) : The phase-space network. Two nodes are connected if they are different by one basic flip (i.e., the flip of a four-spin loop). 

(E) : The flip of a spin loop does not change the height difference, ha — hi, if ha and hi, are both inside or outside the loop. 

(F) : m X n square ice (m = 3, n = 4) with periodic boundary conditions. There are (m — 1) x (n — 1) = 6 types of different 
{Ahx, Ahy). Conflgurations with the same {Ahx,Ahy) form one cluster by adding/removing spheres in the corresponding 
container, while configurations with different {Ah^, Ahy) are disconnected because zero energy flips cannot change the height 
mismatches. Each state has a unique lowest height plaquette labeled in yellow. 
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FIG. 15: (A): 2D circle stacking is combinatorically equivalent to 2D square stacking because each circle or square is supported 
by two circles underneath in a gravity field. (B): Mapping a 2D stack of squares inanax& = 7x9 container to a chain of 
a solid particles and b holes [45]. The dynamics, i.e., the diffusion of particles, is described as a symmetric simple exclusion 
process (SSEP). The number of 2D stack configurations in a container is C"^^ — (a + 6)!/(a!6!), i.e., the number of ways to put 
a particles onto a + b sites. 
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